Stability of self-consistent solutions for the Hubbard model at intermediate 

and strong coupling 



V. Janis 

Institute of Physics, Academy of Sciences of the Czech Republic, 
Na Slovance 2, CZ-18221 Praha 8, Czech Republic 
e-mail: janis(3ifzu. cz 
(February 1, 2008) 

We present a general framework how to investigate stability of solutions within a 
single self-consistent renormalization scheme being a parquet-type extension of the Baym- 
Kadanoff construction of conserving approximations. To obtain a consistent description 
of one- and two-particle quantities, needed for the stability analysis, we impose equations 
of motion on the one- as well on the two-particle Green functions simultaneously and 
introduce approximations in their input, the completely irreducible two-particle vertex. 
Thereby we do not loose singularities caused by multiple two-particle scatterings. We 
find a complete set of stability criteria and show that each instability, singularity in a 
two-particle function, is connected with a symmetry-breaking order parameter, either 
of density type or anomalous. We explicitly study the Hubbard model at intermediate 
coupling and demonstrate that approximations with static vertices get unstable before 
a long-range order or a metal-insulator transition can be reached. We use the parquet 
approximation and turn it to a workable scheme with dynamical vertex corrections. We 
derive a qualitatively new theory with two-particle self-consistence, the complexity of 
which is comparable with FLEX-type approximations. We show that it is the simplest 
consistent and stable theory being able to describe qualitatively correctly quantum critical 
points and the transition from weak to strong coupling in correlated electron systems. 



I. INTRODUCTION 

One of the major problems in the theory correlated electrons is to construct in a systematic and con- 
trolled way a consistent approximation interpolating reliably between weak- and strong-coupling regimes. 
The two extreme limits of weak and strong couplings in the archetypal Hubbard model can be described 
relatively well. The weak-coupling regime is governed by a Hartree-Fock mean field with dynamical fluc- 
tuations covered by Fermi-liquid theory. Extended systems at low temperatures are Pauli paramagnets 
with smeared out local magnetic moments. For bipartite lattices antiferromagnetic long-range order sets 
in at half filling and zero temperature at arbitrarily small interaction. In the strong-coupling regime the 
Hubbard model at half filling maps onto a Heisenberg antiferromagnet with pronounced local magnetic 
moments and the Curie- Weiss law for the staggered magnetic susceptibility, at least at the mean-field 
level. The spectral structure is dominated by separated lower and upper Hubbard bands and the strongly 
correlated system seems to be insulating even in the paramagnetic phase. 

However, it is the intermediate coupling, where the effective Coulomb repulsion is comparable with the 
kinetic energy and hence neither very weak nor very strong, that is of great interest for the theorists as well 
as for the experimentalists. At intermediate coupling, dynamical fluctuations control the low-temperature 
physics of interacting electrons and neither weak-coupling nor atomic-like perturbation theories are ad- 
equate. In this nonperturbative regime a singularity in a two-particle function is approached and we 
expect breakdown of the Fermi-liquid regime and a transition to an ordered state breaking the symmetry 
of the Hamiltonian. Most of the theoretical effort in strongly correlated electron systems concentrates on 
finding an appropriate and reliable description of this fiuctuation-dominated transition region. 

There are only a few general theoretical means with which we could address the weak-to-strong- 
coupling transition. The only exact, Bethe-ansatz solution in d = 1 |^] does not have a weak-coupling 
regime and numerical techniques such as quantum Monte-Carlo simulations or exact diagonalization are 
restricted to relatively high temperatures or small clusters ||^ . They cannot cope efficiently with the very 
different energy scales relevant in the transition regime. Recent progress in the dynamical mean-field 
theory constructed from the limit d oo, where only frequency fluctuations are dynamically relevant, 
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has brought new insights into possible scenarios of the transition from weak to strong couphngs. The 
transition at the mean-field level is connected with a Kondo-like behavior and the Mott-Hubbard metal- 
insulator transition. Its critical behavior has not yet been fully cleared. A widely accepted picture for 
the transition at zero temperature, initially based on the so-called iterated perturbation theory (IPT) ||^ 
and lately supported by analytical and numerical |^] renormalization-group arguments, was, however, 
recently questioned using assessments from perturbation theory and a skeleton expansion j^. Unfortu- 
nately even in this particular limit of high spatial dimensions, no exact solution to the Hubbard model 
has yet been found and our predictions and conclusions are justified by approximate schemes of either 
numerical or analytical character The need for a nonperturbative quantitative scheme being stable 
throughout the whole range of the interaction strength is still extreme. 

We concentrate in this paper exclusively on analytically controllable approximations based on partial 
summations and renormalizations of Feynman diagrams. Only these approaches have the power to give 
the details of the critical behavior of correlated electrons. To decide which of the existing approaches 
are prospective in a reliable description of quantum critical behavior, we investigate in detail stability 
of approximations constructed within renormalized weak-coupling expansions continued to intermediate 
and strong couplings. 

We use the framework of Baym and Kadanoff to formulate general stability criteria that must be 
imposed upon any approximate theory intending to be thermodynamically consistent and conserving. I. e. 
we use an expansion with renormalized one-electron propagators and construct a Luttinger-Ward gener- 
ating functional from which all the physical quantities are derived via functional derivatives. Additively 
we use small external perturbations to test stability of the equilibrium states. 

In this paper we extend the Baym and Kadanoff scheme in two important aspects. First, we generalize 
external perturbations of the equilibrium system from density fluctuations also to anomalous sources not 
conserving charge and/or spin. Only in this way we control all relevant two-particle vertex functions 
that may show divergences and hence lead to instabilities. We use the disturbances as generators of two- 
particle irreducibility and find a complete set of criteria for a local stability of self-consistent solutions 
with a rigorous relation between singularities in two-particle functions, indicating an instability of the 
solution, and symmetry-breaking order parameters. 

Second, since two-particle functions play extremely important role in the stability analysis we have to 
treat the vertex functions and the one-electron propagators on the same footing within a single renormal- 
ization scheme. We must not loose any important relation between various two-particle quantities as well 
as between one- and two-particle functions. This is achieved if instead of approximating the self-energy 
functional we shift diagrammatic approximations to the level of two-particle functions. 

Equations of motion for one- and two-particle Green functions are essential in our construction of 
approximate theories. They are the Schwinger-Dyson and Bethe-Salpeter equations. The former connects 
the one-particle self-energy with the full two-particle vertex function. The latter equations are defined 
separately for each two-particle irreducibility channel and determine finally the full vertex as a functional 
of the completely two-particle irreducible one. Since we expect that the Bethe-Salpeter equations may get 
singular and the full vertex function can have poles at intermediate and strong coupling, we do not attempt 
to destroy the structure of these equations of motion. Instead we introduce diagrammatic approximations 
only in the input to these equations, i. e. in the completely irreducible two-particle vertex. Thereby we do 
not miss possible singularities of the vertex function assumed the relevant physics is essentially determined 
by two-particle scatterings. Such a renormalization of Feynman diagrams, or construction of skeleton 
diagrams, called parquet approach is from the very beginning self-consistent at the level of one- 

as well as two-particle functions. It contains dynamical vertex renormalizations and guarantees that 
divergences in two-particle functions are treated properly. 

The parquet approach or the parquet diagrams are known for long in many-body physics. First 



uickly its way into condensed matter. It 
IJ] , or the formation of the local magnetic 



introduced in nuclear physics [ [12[ , the parquet theory found qi 
was attempted to solve the Kondo problem [|l3| , the x-ray edge 
moment in dilute alloys with the parquet approximation. However, in none of these applications a full 
solution to the parquet equations was found. In the former two approaches an expansion with leading- 
logarithmic divergences of single electron-hole bubbles was used. This construction underestimates the 
role of the genuine poles in the Bethe-Salpeter equations due to multiple two-particle scatterings. It 
actually amounts to a loop expansion neglecting the complex dynamics of the vertex functions. The last 
approach resorted to static approximations and had only a restricted success in the strong-coupling limit. 
It is the intricate structure of the parquet equations that makes nonperturbative solutions inavailable 
and hinders a broader application of the parquet construction. 
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Having a general scheme for deriving theories with renormahzed one- and two-particle functions we 
must reach for approximations to come to quantitative results. The simplest theory fulfilling the equa- 
tions of motion for both one- and two-particle Green functions is the parquet approximation where the 
completely irreducible vertex is just the bare interaction. However, up to now no nonperturbative solution 
to the parquet equations has yet been found. There already exist simplifications of the parquet construc- 
tion (vertex renormalization) jlj] , but because they resort to static vertices they get unstable and fail 
at strong coupling. Here we reduce the intricate complex structure of the full parquet approximation so 
that the vertex functions are dynamically renormalized and allow only for integrable divergences. We 
demand an asymptotic accuracy of the simplification at critical points and keep only maximally diver- 
gent contributions to the vertex function We end up at zero temperature and half filling with 
a theory stable throughout the whole range of the interaction strength the computational complexity of 
which is comparable with the single-channel or FLEX-type approximations. We apply it to a mean-field 
description of the Mott-Hubbard metal insulator transition. 

The layout of the paper is as follows. Section || presents the general scheme for deriving approxima- 
tions with renormalized one- and two-particle propagators allowing us t o formulate stability criteria with 
the relevant two-particle functions. We explicitly show in Sections HI, EV that all approximations sim- 
pler than the parquet diagrams get unstable at intermediate coupling of the Hubbard and single-impurity 
Anderson models. Stability criteria for the parquet approximation are explicitly formulated in Sec. 0. 
The proposed simplification of the parquet approximation with dynamical vertex corrections is presented 
in Sec. VI. We demonstrate stability of the sim plifie d parquet ap proac h in the mean-field theory of the 
Mott-Hubbard metal- insulator transition in Sec. VH. Last Section VIII brings discussion and concluding 
remarks. 



II. EQUATIONS OF MOTION, TWO-PARTICLE GREEN FUNCTIONS AND STABILITY OF 

SOLUTIONS 

We use two generic Hubbard and single-impurity Anderson Hamiltonians with an external magnetic 
field B to model moderate and strong electron correlations. The magnetic field allows us to lift the 
degeneracy in the spin space and to determine the relevant two-particle functions we need to look at if 
we want to prove stability of the theory. We start with equilibrium Hamiltonians 
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The Hubbard Hamiltonian (Rq) is of primary interest for us, since phase transitions with long-range 
order can appear there. The Anderson Hamiltonian (|l^) is considered only to demonstrate the difference 
between a mean-field solution of the Hubbard model (d ~ 00) and the impurity solution and to test 
adequacy of chosen approximate schemes. 

The most direct way to study stability of solutions at thermal equilibrium is to introduce auxiliary 
sources as an external perturbation Only external sources conserving spin and charge leading to 
density-density (nonanomalous) response functions are standardly used to perturb the equilibrium state. 
However, these sources are not complete and do not allow for anomalous responses and quantum phases 
such as superconductivity. To be able to describe all possible instabilities of the Hubbard-like models we 
introduce a generalized external perturbation 
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dld2| [r;ll (1, 2)4(l)c,(2) + (1, 2)c,(l)c,(2) + ^11(1, 2)4(1)4(2)^ 
77^(1, 2)c|(l)c^(2)+ 7,^(1, 2)c|(2)c^(l)l + [e^(l, 2)c^(l)c^ (2) +e^(l, 2)4(2)4(1)1 I (2) 



where labels 1 = (ri,ii), 2 = denote a set of space-time variables characterizing the motion of 

quantum particles, electrons. Here the local real fields 77^ induce density responses conserving charge as 
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well as spin. These fields (e. g. magnetic) are used in standard stability analyses. The new nonconserving 
external sources are complex and either add or remove charge, spin or both from the equilibrium many- 
particle state. The field f]-^ conserves charge but increases spin of the equilibrium state by two elementary 
units. The fields increase charge while increase both charge and the appropriate spin projection. 
The complex conjugate fields lower the respective quantities. 

Fundamental quantity for quantum statistics is a thermodynamic potential as a functional of the 
external perturbation Hext- We need to expand the generating thermodynamic potential only to second 
order in the external perturbation in order to determine (local) stability of the equilibrium state. First 
term of this expansion generates one-particle Green functions. The real fields lead to regular propagators, 
while the complex ones to anomalous Green functions. They vanish at equilibrium unless we are in a 
quantum phase with anomalous Green functions. To decide whether a long-range order may arise or 
whether the equilibrium state is stable we have to evaluate second term in the external perturbation. 
It defines two-particle functions, regular and anomalous. In the unperturbed, equilibrium state without 
anomalous functions only diagonal terms in the external sources do not vanish. If Ha denotes a chosen 



external source perturbing the equilibrium system, i. e. it stands generically for 77II, 77^, 
can write 



and ^ 



we 
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It is interesting to note that second derivatives w.r.t. different external sources introduced in perturbation 
Hext, (^), lead at equilibrium to inequivalent two-particle irreducibility channels, diagrams that cannot 
be made disconnected by cutting two specific lines. More precisely, the external sources generate two- 
particle reducible functions, while their Legendre conjugate "order parameters" the respective irreducible 
functions. The field 77" leads to reducible functions in the interaction ([/) channel and evokes normal 
density responses and susceptibilities. The field rj^ leads to the electron-hole {eh) and to the electron- 
electron (ee) reducible functions, respectively. Second derivative w.r.t. vanishes at equilibrium without 
anomalous functions . 

It is more convenient to introduce a vertex function instead of the full two-particle function G^^^". 
We subtract the free two-particle propagation, if contributes, and stripe off the external legs from the 
remaining two-particle function. The result is a vertex function F independent of the index a [ pT| . The 
full vertex function F can be decomposed equivalently in each two-particle channel into a sum of the 
irreducible and reducible functions A and /C, respectively: 



F = A" 



(4) 



To go on we have to introduce a convention for independent variables labeling the two-particle func- 
tions. We denote X„„i(k,k' ;q) the generic (regular) two-particle function in the momentum represen- 
tation. We use four- momenta, k = (k, ia;„) for the fermionic and q — (q, w„i) for the bosonic variables 
with Matsubara frequencies ujn = (2n -I- l)7r, Vm = 2rmT. The meaning of each variable in the function X 
is manifested in Fig. ^ The incoming arrows denote annihilation of charge. 



ak 



ak + q 




a'k' 



a'k' + q 



FIG. 1. Generic two-particle function with three independent four-momenta and a defined order of incoming 
and outgoing fermions. 



A two-particle irreducible function in a specified channel a can be defined as a functional derivative 
of the self-energy from the a-channel with respect to the renormalized one-electron propagator 
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where the self-energy is a functional of the full one-electron propagator defined from the generating 
Luttinger-Ward functional 

i:"(l,2) = ^M. (5b) 

^ ' ' (5G"(2,1) ^ ' 

Note that even if the anomalous self-energy vanishes at equilibrium the two-particle function is 
nontrivial, since it is conserves both charge and spin. 

Up to now we looked at the equilibrium system from the thermodynamic point of view. We assumed 
we have a Luttinger-Ward generating functional from which we determine all the relevant quantities 
via functional derivatives. However, we do not know how to find the Luttinger-Ward functional. A naive 
way is to choose a particular set of skeleton diagrams directly for the grand potential. A more convenient 
and usual way is to choose diagrams at the one-particle level by fixing the functional S[G']. However, 
in neither of these constructions we have a connection to the fundamental equations of motion for the 
one- and two-particle Green functions. The equations of motion determine the microscopic dynamics 
of quantum systems. If we want to simulate the exact dynamics of the studied system we must try to 
approximate the equations of motion as accurate as possible. 

First equation of motion we have to fulfill is the Schwinger-Dyson equation connecting the one-electron 
self-energy with the vertex function. It reads for Hubbard-like models with a local electron-electron 
interaction: 

^-(^) = Mr E - «^ E k'; q)G„{k + q)G-,{k' + q)G-Ak'). (6) 

^ k' ^ k'q 

The first term on the right-hand side is the Hartree self-energy and the latter contains second and 
higher-order, dynamical corrections. The usual way to deal with the Schwinger-Dyson equation is to 
introduce approximations to the vertex function FfG] on the right-hand side of (||). We thereby get a 
closed functional S[G'] and close the approximation Q. However, at the critical points, where equilibrium 
states get unstable, the full vertex function F is singular. Brute interferences in the Schwinger-Dyson 
equation can make the approximate dynamics inadequate and qualitatively different from the exact one. 
Since we are primarily interested in the two-particle functions and their appropriate description, we 
move approximations from the full vertex function to a lower level, namely to the completely irreducible 
two-particle function. 

Two-particle irreducible vertices determine the full vertex function via Bethe-Salpeter equations. 
These equations explicitly sum two-particle reducible diagrams in each two-particle channel. In the 
notation of Fig. |l| we obtain 

F(fc, k'- q) = A"(fc, k'; q) - [A"GG F] {k, k'; q). (7a) 

The Bethe-Salpeter equations must be completed by definitions connecting the respective irreducible and 
reducible functions A", /C", respectively. Due to inequivalence of different two-particle channels we have 

A" = / + ^ /C"' (7b) 

where I[U;G,T] is the completely irreducible two-particle vertex being a functional of the one-particle 
propagator G and generally also of the vertex function F. When expressed diagrammatically, the com- 
pletely irreducible vertex can be disconnected only by cutting at least three one-electron propagators. 

We used a generic notation for the channel-dependent multiplication of two-particle functions. 
It mixes the variables of two-particle functions in different manners. The three matrix multiplication 
schemes for two-particle quantities in our notation of Fig. ^, representing summations over intermediate 
states, read in the interaction, electron-hole, and electron-electron channels, respectively 



XGG*Y 



{k,k'-q) = ^ E X,a"{k,k";q)G,n{k")G,n{k" + q)Y,„,,{k" ,k';q), (8a) 

~" k'' 
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XGG^Y 



XGGoY 



(fc, k'- q) = ^J2 (fc' fc'; l")GAk + q")G,, {k' + q") 



11 ij , „" 



y.Y,a.{k + q\k' + q"-q-q"), (8b) 



^ (fc, A:'; g) = -L ^ X,,, (fc, k' + q";q- q")G, {k + q - q")G,. (fc' + g") 



-KY,„,{k + q-q",k'-q"). (8c) 
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For convenience we used separate symbols for each multiplication scheme. Note that only the interaction 
(vertical) channel mixes the spin singlet and triplet functions. 

Now all the quantities depend explicitly only on the model input parameters and are functionals of the 
completely irreducible vertex I\U;G,V\to which approximations can be chosen. Two-particle functions 
are determined self-consistently from (0) and (^). The one-electron self-energy is defined from (0). The 
generating Luttinger-Ward functional is determined from the functional differential equation (^). Its 
explicit form for the parquet approximation with I — U is derived in the Appendix. 

Approximate theories break some exact relations and the thermodynamic and dynamical definitions 
of the two-particle functions do not coincide. If we enter an approximation, we have ambiguous def- 
initions for two-particle irreducible functions. They can be defined either from (|^) or from (^. This 
is an unrecoverable feature of all approximate theories which one has to live with. The only thing one 
can demand from reliable approximations is that both the definitions of two-particle functions lead to 
qualitatively the same physical behavior (phase diagram and spectral properties). The vertex functions 
from the equations of motion are responsible for spectral and analytic properties of the approximate 
theory and are used to determine internal stability (consistency) of approximations. On the other hand 
the differential definition of the vertex functions is used for determining thermodynamic behavior and 
stability with respect to external disturbances. They may be either due to external sources or due to 
neglected dynamical vertex renormalizations. 

The advantage of the above general formulation with equations of motion is that we treat the one- 
and two-particle functions consistently on the same footing within one renormalization scheme. This 
is essential for determining stability of approximations. We now rewrite the channel-dependent matrix 
multiplications so as to distinguish manifestly active and conserved variables (parameters) . We introduce 
new notations in which only the active four-momenta are used as variables. We define the following 
symbols 

[XGGf [q]{ak,a'k') = X,,,{k,k';q)G„,ik')G,,{k' + q), (9a) 
[^GG]f,, [q]{k, k') = X,,,{k, k + q,k'~ k)G,{k')G,,{k' + q), (9b) 
[XGG]Z, [q]{k,k')=X,,,{k,k'-q~k-k')G,{k')G,,{q~k'). (9c) 

In (^ the bosonic variable q is always inactive in the multiplication, or conserved during the scatterings 
within the chosen two-particle channel. Each channel, however, has a different conserving variable. It 
makes the parquet algebra of two-particle functions complicated. 

We can, in principle, find eigenvalues and eigenvectors for the matrices, kernels of the Bethe-Salpeter 
equations (0), [A^GG]", and denote them Q". These are complex four-vectors. Only real four-vectors, 
i.e. for zero frequency, both the conserved q" as well as the eigenvectors Q" are important for determining 
the stability conditions. They may formally be written as 

min [A" GG] " [q",0](Q",0) > -1 (10a) 



or equivalently using (7b) 



/ + ^ /C"' I GG 

a' ^a. 



[q",0](Q",0)>-l. (10b) 



This is a complete set of conditions for a local stability of any symmetric solution fulfilling the Bcthe- 
Salpeter equations of motion. The inequalities are a direct consequence of the nonexistence of singularities 
at the real axis in equations (^. Equality in ( p^ indicates that a pole in the respective two-particle 
(reducible) Green function appears and consequently a singularity in a correlation function such as 
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susceptibility. Due to (ffq), the singularity is transferred to the irreducible functions from the other two- 
particle channels and unless the singularity is integrable, the scatterings with singular irreducible kernels 
destroy the transition. We hence can have only integrable singularities in the exact theory. A singularity 
in an approximate two-particle function means on one hand the existence of long-range correlations as a 
precursor of a long-range order, but on the other hand also signals critical vertex corrections stabilizing the 
symmetric phase. These two critical phenomena compete and we must compare the symmetry breaking 
contributions at the critical point with the singular vertex renormalization in approximate theories jl^ . 

The nonexistence of poles at the real axis is a condition for the thermodynamic stability of the sym- 
metric (paramagnetic) phase. However, a consistent solution for S and F must produce analytic functions 
in complex frequencies. This is achieved if the analytic continuations of the two-particle functions from 
Matsubara frequencies have no poles in the complex plane away from the real axis. Hence the proper 
analyticity of the solution is guaranteed if 

[A"GG]"[q",z](Q",Z)^-l (11) 

where lmz,lmZ ^ 0. If these conditions are fulfilled the self-energy determined from (^) is analytic in 
the lower and upper complex half-planes. 

The stability criteria, when applied to approximate solutions, must be used with both two-particle 
functions. The two-particle functions from the equations of motion are used to determine internal con- 
sistency of the approximation and especially ( pT|) is used to prove analyticity of the approximate Green 
functions. The thermodynamically derived vertex functions are used to check stability of the approxi- 
mation with respect to fluctuations caused either by long-range pair correlations or by neglected vertex 
corrections. An "external" instability of an approximation hence need not automatically mean a transi- 
tion to an ordered state. It may happen for oversimplified approximations, what actually is the case in 
low dimensions, that the neglected dynamical vertex renormalization smears out the transition. Hence, 
before we make conclusions on the existence of a symmetry breaking with a long-range order, we must 
be sure the neglected vertex corrections are irrelevant. If not, the approximate LRO becomes unreliable 
and further vertex renormalizations are indispensable. 



III. HARTREE APPROXIMATION 



We investigate in this and the following sections to what degree various approximations produce stable 
equilibrium states at intermediate and strong coupling. Only stable solutions are reliable or suitable for 
the description of the transition from weak to strong couplings. 

We start with the simplest, Hartree approximation. The Hartree approximation is obtained from our 
construction if we put F = in (H). The theory is internally consistent in a trivial way and the only 
nontrivial two-particle function is that obtained from (^ . To determine the external stability we have to 
evaluate single two-particle bubbles 

X„„'{ci,iiyjn) = ■^^Ga(k,iuJn)Gcr'{k + q,i{u;n + Vm)) (12) 

kn 

with the Hartree propagators. They explicitly read for the Hubbard and single-impurity models 

Ga{k, z) = [z + n + aB - Un-a - e(k)]"^ , (13a) 

G^{z) ^[z-Ef+aB- Un^„ - V'^g{z + aB)] (13b) 

where g{z) is the local element of the one-electron propagator of the conduction electrons. We use 
g{z) = 2{z — \J z"^ ~ X^jui^ corresponding to a Bethe lattice in ci = c». 

We evaluate the bubbles for the Hubbard and single impurity models only in the paramagnetic phase 
at half filling, the electron-hole symmetric situation. We start with arbitrary finite dimension. The singlet 
and triplet electron-hole bubbles after summation over Matsubara frequencies are: 

XTx(q, C) = XH(-q, -0 = ^}_. C-2B-Um + e{k)-eik + q) ' ^^^^^ 

Y .^ _ ^ V- / (<tc) - a{B + Um/2)) - f (e(k + q) - a{B + Um/2)) 
X.M. ^^-Jf\. C + <k)-e(k + q) ■ 
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The Hartree solution becomes unstable if 

1 + X,,'(q,0) ==0. (15) 

This happens first for the singlet bubble at a critical field Bcq- It can be understood as a separator 
between the weak and intermediate coupling regimes, since iterations from the Hartree solution get 
unstable beyond the stability limit, Fig. 0. 
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FIG. 2. Instability of the Hartree approximation for the half-filled Hubbard model in an external magnetic 
field for low temperatures in d = 3. The boundary field B^o separates weak and intermediate couplings (WC, 
IC). The energy unit is fixed by the half-bandwidth w = 1. 



The critical field Bco coincides at zero temperature with the boundary to the saturated ferromagnetic 
solution Bg p^ ]. This holds even in the weak-coupling regime with no singularity in the fully saturated 
ferromagnetic state. For if we get below the line of instability of the fully saturated ferromagnet, defined 
hy Bg =11) — U/2, the left-hand side of (j^^ becomes negative and logarithmically diverges with a 
nonvanishing weight p{B -\- Um/2), where p is the density of states and rn < 1. 

We now perform the same analysis for the SIAM in a Bethe lattice and the Hubbard model in d = oo 
dimensions where only frequency variables are relevant and no long-range order can arise. However, the 
strong-coupling phase may get insulating. We choose the following value of the hybridization V = w/^/2 
with the energy unit set by the energy half-bandwidth w = 1. The electron-hole bubble of an impurity 
embedded in a lattice with a finite bandwidth has two contributions, band and impurity ones. They are 
at zero temperature 



dx ^/l — x'^ 1 



r(O)'' = —U / 

^ J -I n U^m^/^ + 1 - a;2 Um/2 + yj{x - 2BY - 1 

U'^m dx Vl-x'^ + Jl-(x-2B)^ , , 

' (16a) 



2 7-1+2B [C/2m2/4+ 1 - a;2] [[/2m2/4+ 1 - (x- 25) 



r(0){, = . (16b) 



^ ' Vm/2 + ^/ ( 25 + ^1 + V^m^ /4 ) -1 

where the local magnetization is 
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Um 



dx 



2^1 + U^m^ /A J-B 



(17) 



Evaluating the instability condition r^(0) + r^(0) = —1 we obtain a critical field beyond which the 
Hartree approximation is no longer stable. It is plotted at zero temperature in Fig. |[ The boundary field 
at which the fully polarized solution (m = 1) breaks down is independent of the interaction strength, 
Bg = w, [^. The Hartree solution becomes exact in the fully spin-polarized state. At i? = 0, the 
instability of the Hartree solution with respect to a spin flip, Uco^ lies below the interaction Up — 2w at 
which the solution breaks into a spin polarized state with m ^ 0. 
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FIG. 3. Instability of the Hartree approximation for the half-filled SIAM in an external magnetic field at zero 
temperature. The dashed line is the local magnetization along the instability line rrico- 



A mean-field solution with only frequencies (i.e. neglecting any long-range order) behaves differently. 
The singlet electron-hole bubble for a Bethe lattice reads 
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with the magnetization defined by an equation 



. , ^ U 
arcsm \ B + —m 
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(18a) 



(18b) 



The boundary field for the fully polarized solution is now a function of the interaction strength alike in 
finite-dimensional models, Fig. ^ The instability of the Hartree approximation outside the fully polarized 
solution is plotted in Fig. ||. Unlike finite-dimensional cases, the weak-coupling regime is robust at zero 
temperature, since we suppressed any possible long-range order. At _B = 0, the instability of the Hartree 
solution appears at Uco = Sttw/S < Up — 'K'w/2. It lies deep below the expected metal-insulator transition 
at Uc ~ iw . It means that the Hartree approximation for effective impurity models gets unstable and 
unreliable before a spin-polarized or insulating solution may appear. We must go to more sophisticated 
approximations iiU > w and we want to approach the expected metal-insulator transition. 
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FIG. 4. Boundary field Bg for the fully polarized solution of the Hubbard model atr = Oonad = oo Bethe 
lattice with local two-particle functions. The boundary field for single impurity is Be = w = 1. 
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FIG. 5. Instability of the Hartree approximation for the half-filled mean-field solution of the Hubbard model 
in an external magnetic field at zero temperature. The upper region (SF) is a fully polarized state (m = 1). 



IV. SECOND ORDER AND SINGLE-CHANNEL APPROXIMATIONS 



An instability in the Hartree approximation in the electron-hole channel, i.e. for the electron-hole 

singlet bubble, indicates either the existence of a symmetry breaking with anomalous functions (in an 
external magnetic field) or an instability with respect to dynamical fluctuations. They can smear out 
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the static Hartree instability and keep the symmetric solution stable. We now check stability of simple 
extensions of the Hartree mean-field theory. 



A. Self-consistent second order 



First step including dynamical fluctuations beyond Hartree is second-order perturbation theory. We 
obtain this approximation if we replace the full vertex function T in by the bare coupling constant 
U. Since the right-hand side of (^) contains no infinite series, the only stability criterion for this approx- 
imation is that with (|5a|). Alike the Hartree approximation we expect that the instability arises first in 
the electron-hole channel with spin singlet bubbles 

eh 

[q]ik, k') = U{1 - X,_M)Ga{k')G^„{k' + q) (19) 

cr — cr 

where Xaai{q) was defined in (]T^). We have to find the lowest real eigenvalue of the above expression 
seen as a matrix in k and k' . Since the matrix does not explicitly depend on k, the lowest eigenvalue is 
just a sum over the other variable k' . We then have a stability criterion 

X,_,(q,0)(l-X,_,(q,0)) > -1. (20) 

The linear term on the left-hand side origins from the Hartree approximation. The other term 
from second-order theory formally enhances the tendency to instability of the Hartree solution, since 
Xo--[T(q, 0) < 0. However, the self-consistence renormalizing the one-electron propagator may push the 
instability to slightly higher values of the interaction strength. 



GO 



B. Iterated perturbation theory 

Iterated perturbation theory is based on the impurity second-order perturbation theory |24p[|. It 
seems to be up to now the only simple diagrammatic theory capable to describe both the metallic (Fermi 
liquid) as well as the insulating (atomic limit) phases. This mean-field theory is partially self-consistent 
with frequencies as dynamical variables. A correction to the Hartree self-energy ASo- within IPT reads 

AE<^(ia;„) = -^'^Qa{i^n + Wm)X^a^a{iVm) (21) 
m 

where 

Q^i^)= ,^rfl , , (22) 

1 + Ga{z)T.a{z) 



and A" is a two-particle bubble defined as in (|12|) but with the frequency-dependent (local) propagator 
Q . We define a new two-particle quantity 

t/ 1 1 

X / dtp{e) ^- ^ %— ^ (23) 

needed for the stability criterion of IPT solutions. We easily find that the IPT solution is stable w.r.t 
spin flips if 

^.-.(0) - X^-M [(^-,-.(0)) - X,-AQ)] > -1. (24) 

We see that IPT is generally more stable (to higher values of the coupling constant) than the fully self- 
consistent perturbation theory from the preceding subsection. It is due to subtraction in the parentheses. 
This difference stems from a specific dependence of the one-particle propagator Q on the self-energy. The 
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closer to the impurity case we are (narrow energy band) the smaller the dynamical contribution to the 
instability of the IPX solution is. In the extreme case of the impurity model, there is no contribution 
to the instability from second order, i.e. beyond Hartree. Because of the Hartree term there is always a 
critical interaction at which the IPT solution gets unstable in both the single impurity and the mean-field. 
This instability, as we already know from the Hartree solution, lies well below the expected Mott-Hubbard 
metal-insulator transition. It is hence inappropriate to extrapolate a metallic IPT solution to the critical 
region near the Mott-Hubbard transition and beyond. It seems to be no way of matching consistently 
the strong-coupling (insulating) and the weak-coupling (metallic) IPT solutions. 

IPT from mean-field models amounts to a non-self-consistent expansion on finite dimensional lattices. 
A non-self-consistent expansion is derived within the general scheme of Section || if the renormalized one- 
electron propagator is replaced by G^°-'^^ = G^^ + T,[G\ where the self-energy is treated as a functional 
of the full one-electron propagator. All orders of the non-self-consistent expansion get unstable at the 
same point where the Hartree approximation does. A non-self-consistent expansion does not improve 
upon stability of the Hartree solution. 



C. Single-channel approximations 

Dynamical approximations cannot push the instability of the Hartree solution to much higher inter- 
action strengths unless they contain infinite-many two-particle bubbles. Simplest theories with geometric 
series of two-particle bubbles are single-channel approximations or their extension, the so-called FLEX 
approximation |^ . I. e. we explicitly sum multiple two-particle scatterings on the bare interaction within 
a single two-particle channel. We derive these approximations within the general scheme of Section ^ if 
we choose A" = [/, A" = 0, a' ^ a Since we expect that the instability arises in the electron- hole 
channel, we explicitly consider only a renormalized RPA (electron-hole channel). A correction to the 
Hartree self-energy for this approximation is 

AE,,.) = -iLj^c_.» + „^M_^ (25) 

The full vertex function F from the Schwinger-Dyson equation (||) already contains a geometric series of 
two-particle bubbles, so we have two nontrivial stability criteria. The one derived from the equation of 
motion in this approximation reads 

^a-a(q,0)>-l. (26) 
This is the stability condition from the Hartree approximation with a renormalized one-electron prop- 



agator with the self-energy from (E9). Condition (26) is fulfilled for weakly and moderately correlated 



systems significantly beyond the instability of the Hartree solution. In the impurity or mean-field cases 
it is even satisfied up to infinite interaction strength H^] . A solution fulfilling (|2^) is internally consistent 
and possesses the necessary analytic properties of the self-energy. It, however, does not guarantee stabil- 
ity of the approximation with respect to external perturbations and to fluctuations beyond this solution. 
For that purpose we use definition (^^ for the irreducible functions. With the generic notation (^) we 
obtain a new vertex function in the electron-hole channel 



5G 



-\ eh 

[q]{k,k') 



u 



l + X^_^(fc'-fc) 

f/2 ^ G^^ik + k' + q- k") G„{k") 



u \ ^ 

'^STV ^ r+ X^_^(fc' + q- k") 1 + X„^„{k' - k") ■ 



(27) 



This vertex function enters the irreducible functions in the other channels, the interaction and the electron- 
electron ones. We identify in each channel the conserving bosonic variable according to rules (^) and 
obtain a matrix in k, k' that is to be diagonalizcd. Unlike second-order we cannot find the lowest real 
eigenvalue in closed form. We can at least qualitatively assess the lowest eigenvalue in the asymptotic 
limit of condensation of electron-hole pairs, where the electron-hole two-particle vertex function diverges. 
First, the electron-electron channel, at least at half filling, is irrelevant and we resort to the electron-hole 
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one. Next, we neglect the convolution of the singular vertices with regular electron- hole bubbles assuming 
that the integrals smear out the singularities. We hence keep only the first term on the right-hand side 
of ( |27| ) as a leading divergent contribution in the asymptotic limit of the critical point. Even then we 
still have a full fc, k' matrix to diagonalize. We further resort to zero temperature with no difference 
between fermionic and bosonic variables and assume that the lowest eigenvalue can be approximated as 
in second-order theory. I. e. we treat the matrix as if it were fc-independent. This assumption does not 
hold strictly but only approximately if the influence of the states away from the Fermi surface is not too 
big. Within this approximation we obtain a new approximate thermodynamic stability criterion in the 
interaction channel 

U ^ G-.((ko, 0) + gOG_.((ko, 0) + q' + (q, 0)) 

where ko is a suitable vector from the Fermi surface minimizing the expression on the left-hand side. 

We now have two criteria that both must be fulfilled in order to have a stable phase with correct 
analytic properties of the Green functions. The former criterion, (p6|), from the equation of motion 
guarantees internal consistency and analyticity of the approximation. The latter one, (p8|), when fulfilled, 
implies stability of the approximation with respect to fluctuations beyond the renormalized RPA. When 
broken, scatterings in the neglected interaction channel become relevant and must self-consistently be 
interwoven with the electron-hole ones. Simple additions of channels as in the FLEX approximation do 
not introduce a self-consistence in the two-particle functions and do not improve upon the stability of the 
single-channel approximations. In low dimensions {d = 2, 1), the left-hand side of (^8|) may even diverge 
and the single-channel approximations become qualitatively wrong, since the full solution can produce 
only integrablc singularities. Hence FLEX- type approximations are unsuitable for the description of the 
weak-to-strong-coupling transition. 



V. PARQUET APPROXIMATION 



The approximations we have hitherto dealt with were either derived by simplifying the full two- 
particle vertex function F in the Schwinger-Dyson equation or some of the Bethe-Salpeter equations were 
disregarded. The various two-particle functions were hence treated inequivalently and some divergences 
went lost. The simplest approximation consistent with all the Bethe-Salpeter equations (0) is the parquet 
approximation with the unrenormalized completely irreducible vertex I[U;G,T] = U. We formulate 
explicitly the stability criteria for this approximation and show how the singularities from the Bethe- 
Salpeter equations are controlled. 

To control the singularities that may arise in two-particle functions we diagonalize the Bethe-Salpeter 
convolutive equations (Q). We, however, can diagonalize the Bethe-Salpeter equations only separately in 
each channel and cannot express the solution in closed algebraic form. Diagonalization in the chosen two- 
particle channel means diagonalization of the appropriate matrix from (^. To formalize this procedure 
we represent the one- and two-particle functions as operators on a Hilbert space of two-particle (bare) 
states. Such a Hilbert space can be spanned on the basis vectors |fc, /c' > parameterized by two one-particle 
four-momenta. The two-particle vertex functions can be represented as 

X^^, = S{k + k' -k -k')X^^,{k,k'-k,k')\kk'){k'k\. (29) 

We now switch from this channel-independent representation to channel-dependent diagonal representa- 
tions. We use q for the conserved four-momenta and Q for the eigenvalues of the corresponding matrices 
in the two-particle channels from (^. 

We need diagonal representations for the a-reducible two-particle functions 

Z" = Y^JC^{q,Q)\Qq,a){a,qQ\. (30) 

q,Q 

Note that the eigenvectors \Qq, a > are expressed as a combination of the vectors \k,k' > where Q is a 
combination of fc's. The vector q fixes either difference (electron-hole and interaction channels) or sum 
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(electron-electron channel) of k',k. All \k,k' > and \Qq,a > form complete orthogonal sets. For the 
sake of simplicity we suppressed normalization factors in the sums in ( p9| ) and (pO|). Representation ( |30| ) 
makes the Bethe-Salpeter equations algebraic and we can write their formal solutions as 



{a,qQ\ X« \Qq,a) 



l + {a,qQ\ \Qq,a) 

where X and Y are two-particle operators. They have the following explicit representation 

{a,qQ\ X" \Qq,a) = {a,qQ\UG(^^'>U\Qq,a) + ^"'(^''Q') 

X Ya,qQ\UG^^'^\Q'q',a'){a',q'Q'\Qq,a) + {a,qQ\Q'q\a'){a' ,q'Q'\G'^''^U\Qq,a) 

+ E E /c"'(g^Q0^"''(9'^o'')(«,'?QlQV,a'>(«^9VlG(2^lQ'v^«'': 

a'a"^a q' q" ,Q' Q" 

x{a",q"Q"\Qq,a), 

and 

(a, gQI IQg, a) = (a, gQIf/G'^) |Qg, a) + ^ ^ /C"' (g', Q') 

X (a, <7Q|C/G(') IQ'q', a') (a', q'Q'\G<-^^ \Qq, a) . 
We used abbreviations for the free two-particle propagation 

Gi^j, =^G.(fc)G.,(fc')|A:fc')(fc'fc| 

kk' 

and for the operator of the bare vertex 

U^^,^US„,,-^ S{k + k'-k-k') \kk'){k'k\. 



(31a) 



(31b) 



(31c) 



(32a) 



(32b) 



The parquet approximation has become a set of closed equations for IC°'{q,Q) for the given one-electron 
propagators. The most difficult task in solving parquet equations (pi]), ( |3^ ) is the diagonalization in 
the separate channels, i. e. finding the bases \Qq,a > represented in the bare states \kk' >. They 
are determined from (^ an d (|7b| ) so that the parquet equations remain implicit and must be iterated. 
However, parquet equations (|3l|) allow us to control the behavior of the denominators, i. e. the appearance 
of singularities. 

To close the parquet approximation we have to rewrite the Schwinger-Dyson equation with the diagonal 
form of the vertex function. It reads 



{k'k\UG^^lu\kk') + E E ^-^^(9' 

a q.Q 

x{k'k\Qq,a){a,qQ\UG'^2jkk') 



(33) 



The parquet approximation is completely determined by numbers T,a-{k), /C"^/(g, Q) and < a,qQ\kk' >. 
The internal stability criteria for the parquet approximation in the diagonal representation are 



l + (a,(q,0)(Q,0)| Y" |(Q,0)(q,0),a) >0. 



(34) 



We see that all the two-particle functions can contain only integrable singularities and that if a singularity 
arises in one channel it cannot cause a havoc in the other channels, since the integral kernels are integrable 
functions. 
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There are also stability criteria for the parquet approximation with respect to external perturbations 
beyond this approximate equilibrium state. Since possible instabilities must be integrable we do not 
expect a qualitatively different behavior. A qualitative change could be caused by new bound states with 
more than two particles that could not be decomposed into simpler two-particle states. It is physically 
plausible to assume nonexistence of such states, at least all the existing approaches do that. 

It is a tremendous task to solve the complete parquet approximation. It can be completed only in 
a very limited special situations p^ . However, the general scheme how to solve the parquet equations 
helps us to identify the principal features any simplification should not loose. It is a two-particle self- 
consistence allowing only for integrable singularities and an access to the diagonal or algebraic form of 
the solution of the Bethe-Salpeter equations enabling to control possible singularities. 

VI. SIMPLIFICATIONS IN THE PARQUET APPROACH 

Parquet diagrams differ significantly from the simpler single-channel approximations with multiple 
two-particle scatterings. The difference manifests itself in the effective interaction in the two-particle 
scattering processes. The bare interaction remains unrenormalized in the single-channel approximations. 
In the parquet equations, it is replaced by renormalized vertex functions representing a dynamically 
screened inter-particle potential. A qualitative difference between a solution to the full parquet equations 
and to the single-channel approximations becomes evident at the two-particle criticality. The vertex 
functions from one or more channels get critical and sharply peaked around the Fermi energy. The 
exchange between the scattered particles becomes strongly energy dependent and cannot be approximated 
by a static effective interaction as in the single-channel theories. 

The Bethe-Salpeter nonlinear integral equations are extremely complicated with an intricate analytic 
structure of the solution. Each two-particle function entering the parquet equations contains initially 
three independent (complex) variables. Neither of them can simply be neglected, since the variables are 
interconnected due to a mixing of different inequivalent channels. Even the diagonal form of the parquet 
equations with only two relevant variables remains implicit and does not allow for a solution in closed 
form. That is why a nonperturbative solution to the parquet equations has not yet been found. From this 
reason the parquet diagrams are relatively little used in condensed matter although being in existence 
more than thirty years. 

The only hope that the parquet approximation may become useful for quantitative assessments of 
the effects of strong electron correlations is that one succeeds in simplifying the full parquet algebra to 
an analytically manageable form. As already discussed, the parquet approach is expected to produce 
qualitatively new results at intermediate and strong coupling with critical two-particle functions where 
weak-coupling approximations get unstable. In simplifying the parquet equations we will demand asymp- 
totic accuracy at the two-particle criticality, particularly for condensation of electron-hole pairs leading 
at half filling to a metal-insulator transition. We hence retain only the diagrams potentially making the 
effective mass of the electrons divergent at the metal-insulator transition. 

A. Two- Variable Two-Particle Functions 

It was recently argued that at least at half filling only two of the three topologically inequivalent 
channels are relevant in the critical region of the metal-insulator transition. Only multiple scatterings in 
the electron- hole and the interaction channels contain divergent diagrams at the critical point |Q . This 
simplification to two channels itself, called dipole approximation, does not reduce the number of variables 
in the two-particle functions. It only simplifies the algebra of the parquet approximation. To lower the 
number of the relevant variables we use the fact that the singularities in the two-particle functions in the 
parquet approximation must be integrable. It means that when integrating over the variable in which 
the singularity arises we obtain a finite result bounded by a constant of order one. Our approximation 
consists in neglecting all finite (nondiverging) contributions at the critical point. This must be done at 
the level of the diverging two-particle functions. 

As a first step we reduce the number of the relevant variables in the two-particle functions to two. This 
is achieved when we neglect mixing of the reducible functions from different channels. We replace the full 
two-particle function F on the right-hand side of the Bethe-Salpeter equation in the a-channcl by + 



15 



i. e. by a sum of the completely irreducible vertex and the reducible function in that channel. Such a 
replacement does not change the denominator of the diagonal representation of the solution to the Bethe- 
Salpeter equations and hence the universal critical behavior of the function K."' remains unaltered. The 
bare interaction is a constant and consequently the solution does not depend on the outgoing (incoming) 
variables. Further on, the triplet functions can be eliminated from this approximation, since Ag-.o- gets 
irrelevant. If we denote K"^ _^ = U + ^^'^-cr ^^^"^ A^ _^ = C/ + K,^ the parquet equations reduce to 

-JmT. l')GAk + q')G-,{k + q' + q)K,-Ak + q) , (35a) 

«' 

A:j^_,(fc; q) = U+ q')G-Ak + q')G^,{k + q' + q) 

q',q" 

xA"_,Jk + q'; q" - q')GAk + q")GAk + q" + q)A'^,-Ak + <z"; <z) • (35b) 

Because of the broken symmetry between the incoming and outgoing variables in the two-particle 
functions we must be careful in simplifying the equation for the self-energy. To be consistent we use a 
symmetrized formula with the bare interaction at the incoming and outgoing momenta. We obtain 

AS.(fc) - E + q)G^,{k')G^,{k' + q) K,_,(fc + q; -q) 

k' ,q 

+At^Jk'; q) + Al^_,{k + q;k'-k)+ A%Jk'; k ^ k') - 2U] . (36) 

The reduced two-channel parquet approximation with two-variable two-particle functions reproduces 
in leading order the critical behavior of the parquet approximation (with two channels). The universal 
quantities derived from the divergent functions of the parquet equations are reproduced by equations (^5|)- 
(|3^). Only nonuniversal quantities such as the critical interaction do not come out in this simplification 
as in the full parquet theory. 

We achieved a simplification of the full parquet algebra without loosing leading critical behavior of the 
solution, but the resulting equations retain a nonlinear convolutive character. It is still too elaborate to 
reach nonperturbative solutions or to diagonalize them. We must find a further reduction of complexity 
of the parquet equations to make them useful for an affordable quantitative analysis. 



B. One- Variable Two-Particle Functions 



It is the bosonic variable q in the two-particle functions A'', A" that is relevant at the critical point 
and in which the singularity aris es. T h e fer mionic variable k in both functions "labels" the eigenvalues of 



the integral kernel in equations (35a), (35b). The minimal eigenvalue governs the critical behavior at low 
temperatures. If we are interested only in the asymptotic critical behavior we can approximate the actual 
minimal eigenvalue by a value at the lowest-lying fermionic four-momentum. By such an approximation 
we neglect mixing of the fermionic four-momenta in the Bethe-Salpeter equations and fc-dependence in 
the two-particle functions. This approach corresponds to a low-energy expansion used by Hamann 
by assessing the strong-coupling limit of the rcnormalizcd RPA of Suhl in the single-impurity problem. 
We assume that the approximate minimal eigenvalue dominates and qualitatively determines the relevant 
physics of the critical point. 

The suggested simplification is possible only at zero temperature, where the difference between the 
fermionic and bosonic Matsubara frequencies vanishes. We hence put T = 0, which is the most interesting 
case for the metal-insulator transition. Using the above ansatz, neglecting fc-dependence in A'', A", and 
putting fc — (ko,0) in the one-electron propagators we reduce the parquet equations to a manageable 
form. We introduce new functions 

= ^ E A^,_,(<z')G.((ko, 0) + -z')G-.((ko, 0) + g' + q) (37a) 
^-('?) = ^E^---(9')^-((ko,0) + -7')Ga((ko,0) + g' + g) , (37b) 
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with the aid of which we obtain for the two-particle vertex functions 



The introduced fermionic momentum (ko,0) from the Fermi surface is an adjustable parameter chosen 
to maximize the tendency to instability of this approximation. 

The above equations have the complexity of the single-channel approximations, but contain renormal- 
ized vertex functions in the two-particle scatterings. It is always the renormalized vertex function that 
enters the denominator of equations (|3^) via the functions V and /C and the critical behavior is properly 
renormalized. 

It is easy to derive the corresponding equation for the self-energy. It reads 

AS.(fc) = E + ^)X-a.-.{q) (A^.,(<z) + A^,_,(-g) - U) 

q 

+G^Ak + q)X,,-.{q) {Kl-M + I-.A-l) - U)] (39) 
where Xa-^a' are two-particle bubbles with the renormalized one-electron propagators. 



Equations (|37|)-(39) replace the original parquet approximation. They fully determine the generating 
two- and one-particle functions. All thermodynamic and spectral quantities can be calculated from them. 
Since we used a number of approximate steps in the derivation of the reduced parquet algebra, any new 
quantity must first be formulated within the full parquet approximation. The Luttinger-Ward functional 
generating all thermodynamic quantities in the parquet approximation and simplifications thereof is 
constructed in the Appendix. Simplifications are introduced in the exact formulas via the generating 
one- and two-particle functions. 

Having a simplified and manageable set of parquet-type equations we can finally formulate internal 
stability criteria for them. They are derived from the general ones and resemble the stability criteria 
from the single-channel approximations. They read 

1 + /C,(q,0) > 0, l-fr,(q,0)>0. (40) 

Unlike the single-channel approximations, stability criteria ( ^0|) contain renormalized dynamical two- 
particle vertex functions that are determined self-consistently from the simplified parquet equations. The 
functions To- and JCa depend implicitly on the momentum kg from the Fermi surface emerging due to 
consistency of our simplification. The vector ko is chosen so that to minimize the left-hand side of (|40| ) 
The vector then determines the character of the symmetry breaking (s vs. d wave etc.) first appearing 
in the system. 



VII. DYNAMICAL MEAN-FIELD APPROXIMATION AND METAL-INSULATOR 

TRANSITION 

We apply the dipole (simplified parquet) approximation to effective impurity models (single impurity 
or dynamical mean-field) at intermediate and strong coupling. The fluctuations in space do not contribute 
to the dynamics and the four- momenta from the general formulation of the parquet approximation collapse 
to frequencies. This simplification enables one to perform explicitly analytic continuation from the 
Matsubara to the real frequencies using contour integrals. Realizing that the reduced equations hold for 
r and introducing C, z for the bosonic, fermionic complex frequencies, respectively, we can write for 
the two-particle functions 



icAC) = -ul ^ 



G,(c^ + C)Im 



GaiuJ-C) 



i + r_,(t^+) i + r_,(c^-c) 



lmGaiuj+) 



(41a) 
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r.(C) 



-u 



duj 



G_^(w + C)Ini 



GALU+) 



1 - /C<j(cj+)/C_o-(cJ+) 



1 - /Co- (o) - c)/c_o-(t^ - C) 



ImG'_cr(w+) 



(41b) 



and analogously for the one-particle self-energy correction to the Hartree term 
AEo(z)= r —{GA^ + z)lm[X^,.^Auj+){A^_^Juj+)^U)] 

J — OO ^ 

+X_,,_o(c^ - z) {A''_,Ju;+) - U) ImGo(cc;+) + G^^uj + z)Im [X,,_o(w+)A^^_,(cc;+)] 
+Xo,-<x(t^ - z)kl^^Auj - z)ImG_o(w+)} 



(42) 



where we denoted uj^ = lj + iQ^ . 

Equations (|4^) and (|4|) can be solved numerically. We use the one-electron propagators from the 
Bethe lattice in d = oo with a finite half-bandwidth w = 1. We as well put i? = 0, since the Mott-Hubbard 
metal-insulator transition is expected there and both channels should show the same divergence at this 
critical point. 

We have three levels of self-consistence in the parquet equations. We have to determine self-consistently 
the Hartree parameters (spin-dependent particle densities ), th en the dynamical self-energy correction, and 
finally the two-particle vertex functions (cf. Appendix, (A. 7)). It appears that we have to perform the 
self-consistences sequentially from top to bottom. First the Hartree must be self-consistently determined, 
second the self-energy correction, and finally the vertex functions. Only in this order we can secure 
convergence of the iteration scheme and analytic properties of the converged and intermediate solutions. 

First step in our iteration scheme for the two-particle functions reduces to self-consistent second order. 
It becomes unstable at around U ~ w and gets spin polarized at intermediate coupling. Iterations to the 
FLEX approximation with two (electron-hole and interaction) channels converge very quickly at weak 
and intermediate coupling. The converged solution never gets spin-polarized. However, starting with 
U ~ w the solution becomes unstable with respect to external perturbations, spin- flips, and condition 
( p^ ) is violated. Fig. ^ shows the real part of the vertex function F for the FLEX with the first next 
iteration towards the parquet solution. A stable solution demands that both functions at = lie above 
the value —1. 
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FIG. 6. Vertex function F calculated within the FLEX (solid line) and first iteration beyond FLEX (dashed 
line). Their values at zero frequency decide about stability of the approximation (Rer(O) > —1). 
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The iterations for the reduced dipole approximation converge at weak and intermediate couphng for 
U < 2w. We experienced problems with convergence for higher values of the interaction strength [ p9[ . 
Due to logarithmic singularities we had insufficient resolution close to the Fermi energy. We used a linear 
distribution of the frequency mesh points with Auj = 0.005w, since it is convenient for evaluating the 
convolutions. We need all the energy scales to reach a solution fulfilling necessary sum rules, but as the 
metal-insulator transition is approached, a broad spectrum of energy scales is to be covered. A more 
sophisticated method for dealing with the different energy scales in equations (|4^), such as numerical 
renormalization group, has to be applied if we want to continue the approximation closer to the metal- 
insulator transition. 

Although we were unable to iterate the solutions for interactions approaching the metal-insulator 
transition, we nevertheless can demonstrate a qualitative difference in the behavior of the FLEX and the 
parquet approximations. Apart from the fact that the FLEX approximation becomes unstable, the main 
qualitative difference is evident for small frequencies. Fig. ^ shows the difference in the frequency behavior 
of the real part of the vertex function T at U — 2w. A tendency to forming of a narrow quasiparticle peak 
around the Fermi energy is evident in the parquet approximation. The parquet approximation is much 
closer to the instability point r(0) = — 1 than the FLEX solution. However, due to the self-consistence 
in the parquet approach, the solution remains stable. 




FIG. 7. Comparison of the two-particle vertex calculated at f/ = 2w within the FLEX and the simplified 
parquet approximations. The inset shows the detailed structure around the Fermi energy. 

Fig. H and Fig. |^ show the same for the imaginary part of the one-electron propagator and the self- 
energy. A tendency to the expected insulating solution and splitting of the quasiparticle central band 
at low frequencies can again be clearly recognized. However, the low-energy dynamics in the parquet 
approximation is different from the unstable IPT solution |7j. The central quasi-particle peak is much 
more narrow, however much less separated from the excited states. No pronounced quasigap is yet 
present. The high frequency behavior of the density of states should be dominated by satellite peaks. 
Emerging of shoulders for the expected lower and upper Hubbard bands can be seen in the parquet 
solution (but not in FLEX) . Also the satellite bands (the high-energy region) are less pronounced than in 
the IPT and numerical renormalization group jsj or Monte Carlo simulations ^,0. Hence the separation 
of low- and high-energy scales seems to go much slower in the diagrammatic approaches than in numerical 
simulations. It can be explained by the fact that the spin-symmetric atomic limit with a discrete energy 
spectrum is not easily reproducible in the parquet construction with a continuous spectrum of energies. 
It is the half-bandwidth w that fixes the energy scale and alike the Bethe-ansatz solution it cannot simply 
be limited to zero. No tendency to a divergence in the self-energy is observed, as seen in the numerical 
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renormalization group ||5|]. It is impossible to produce a singular self-energy in the metallic phase in the 
Schwinger-Dyson equation with integrable two-particle functions. 



o 




FIG. 8. Comparison of the imaginary part of the one-electron propagator from the FLEX and the simplified 
parquet approximations. The inset shows a tendency towards the formation of a central quasiparticle peak in the 
parquet approximation. 
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FIG. 9. Comparison of the imaginary part of the self-energy with the layout as in Fig. ^. 

The parquet approximation has two vertex functions that are to show the same divergence. However, 
their dynamical behavior beyond the Fermi energy is different as shown in Fig. ^ It is the irreducible 
vertex function from the interaction channel JC that has a stronger tendency towards instability. The 
values of the two vertex functions at the Fermi level coincide within the numerical resolution and lie very 
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close to —1, within the precision of 10 




The critical point of the metal-insulator transition is not directly accessible with the present numerical 
solution of the simplified parquet approximation. However, the low frequency behavior of the divergent 
functions can be assessed in leading order analytically. At the critical point of the metal-insulator 
transition both vertex functions, U /{I -f r(cj)) and C//(l + K.{lo)), diverge at the Fermi energy w = 0. 
Consequently, the derivatives of the real part of the self-energy (effective mass) and the derivatives of the 
imaginary parts of the vertex functions diverge. It is straightforward to obtain in leading order 

TT l+l_cr(Oj TT 1 - /C_cr(0j/Cc(0j 



Re Y.'M - - - 

TT 



Im G-ct(iO )- -I- Im Ga\ivi 



i + r,(o) ^i-/c_,(o)/c,(o) 



(43b) 



where the small (Kondo) energy scale is A = w{\ + r(0)) or equivalently A — w{l + /C(0)). We 
see from (^) that the divergence of the effective mass is bound to a divergence in the electron-hole 
dominated irreducible vertex functions. It is important to stress that a singularity in one channel causes 
a nonanalyticity in the other channel. Hence suppressing any of the channels with electron-hole loops 
in the effective impurity problems means inappropriate treatment of the metal-insulator (Kondo) critical 
region. 

If the Schwinger-Dyson equation of motion is fulfilled the effective mass cannot diverge, or a sharp 
metal-insulator transition cannot appear without a singularity in the electron-hole (local) vertex functions. 
Fig. |l^ shows for [/ = 2it; the three functions the slopes of which should be infinite at the transition point. 
The least sharp slope displays the vertex function F as can be expected from the more enhanced tendency 
to instability in K. than in F. 

The spin-symmetric situation is the most difficult to analyze because of the degeneracy in the spin 
space (singlet and triplet scatterings are indistinguishable). If we switch on an external magnetic field 
we remove the spin symmetry and the density of states will no longer be pinned at the Fermi level. The 
value oi K,cr{Q)lC-a{Q) decreases and we push the solution away from the criticality of the metal-insulator 
transition. Going from weak to strong coupling at a fixed field, we gradually increase the absolute value 
of Fo-(O) as long as we reach a critical point (for sufficiently small fields) Fo-(O) = — 1 at the boundary 
to the fully spin-polarized (Hartree) state Us- The value of /Co-(0)/C_ct(0) first increases with interaction. 
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reaches a maximum at intermediate coupling and then, due to the Hartree band-sphtting terms in the 
self-energy, decreases to zero at the boundary to the fully spin-polarized state. There is no indication for 
a Mott-Hubbard metal-insulator transition in an external magnetic field. An extrapolation from finite 
fields to i3 = suggests a hypothesis that the metal-insulator critical point is shifted to infinity as in 
the single impurity case. However, neither finite fields do allow for a stable numerical solution with 
linearly distributed mesh points at strong coupling to support the hypothesis quantitatively. A detailed 
quantitative analysis of the Hubbard model in an external magnetic field will be presented elsewhere. 
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FIG. 11. Three functions with divergent slopes at the metal-insulator transition from (0) at U — 2w. 




VIII. DISCUSSION AND CONCLUSIONS 



We demonstrated in this paper the importance of a proper formulation of approximate theories in 
order to determine the range of their applicability. We showed that a careful stability analysis is needed in 
intermediate or strong-coupling regimes due to instabilities, poles that may arise in two-particle functions. 
We must have a reliable description of two-particle functions to be able to decide about stability of 
solutions. We suggested that the parquet approach formulated here as an extension of the Baym-Kadanoff 
construction with properly chosen external disturbances offers the desirable framework for a consistence 
and stability analysis of approximate solutions. Under a physically plausible argument that singularities 
in two-particle functions are due to poles in the Bethe-Salpeter equations, the parquet approach leads to 
a qualitatively correct asymptotic critical behavior at transition points, ft introduces a new two-particle 
self-consistence causing that two-particle functions can contain only integrable divergences as in the exact 
theory. 

We found a complete set of criteria for local instabilities from the poles in the Bethe-Salpeter equations. 
We attached to each singularity in two-particle functions an order parameter breaking a symmetry of the 
underlying Hamiltonian. The parameters are new self-energies and are either of density type, such as 
magnetization or have anomalous character and do not conserve either spin (transverse magnetic order) or 
charge (superconductivity). Since any approximate self-consistent theory breaks some relations between 
various Green functions, we always have two sets of two-particle functions and stability criteria. The two- 
particle functions determined from the Schwinger-Dyson equation of motion are used for internal stability 
and consistency of the theory and determine whether spectral functions have the desired analytical 
properties. The two-particle functions defined thermodynamically from the Luttinger-Ward functional via 
functional derivatives are used to determine external stability w.r.t. external disturbances or dynamical 
fluctuations. The two different two-particle functions must produce qualitatively the same phase diagram 



22 



and spectral properties in a trustworthy theory. 

Approximations in the Baym scheme are introduced for the full two-particle vertex function in the 
Schwinger-Dyson equation. Since the vertex function may get singular, approximations to the Schwinger- 
Dyson equation can make the resulting theory unreliable at two-particle criticality. The presented parquet 
approach does not touch the Schwinger-Dyson and Bethe-Salpeter equations of motion and introduces 
approximations only in the input to these equations, namely in the completely irreducible two-particle ver- 
tex. This is a new qualitative change in describing critical behavior, since unlike the full vertex function, 
the completely irreducible vertex is not expected to cause new singularities to appear. Diagrammatically, 
it can only be broken into two separate parts by cutting at least three one-particle propagators. It reduces 
in the parquet approximation to the bare interaction. A lowest-order correction to it is plotted in Fig. [l^ . 
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FIG. 12. A lowest-order contribution to the completely two-particle irreducible vertex beyond the parquet 
approximation. The bare interaction is replaced by the full two-particle vertex from the parquet solution. Double 
primed variables are summation indices. 



Corrections to the completely irreducible vertex beyond the parquet approximation consist of integrals 
over potentially singular functions and do not generate new divergences. These corrections may, however, 
change the critical exponents ]3l| ]. However, an expansion beyond the parquet approximation must 
not be interchanged with a weak-coupling loop expansion. The parquet approximation can be viewed 
upon as a sophisticated "mean-field" theory with dynamical vertex corrections for strongly correlated 
systems. Strongly correlated electrons do not allow for simpler static mean-field theories, since two- 
particle criticality may realistically be described only if we do not loose integrability of divergences and 
stability in all two-particle channels. All simpler approaches loose this essential property of the exact 
theory. 

The most difficult problem with the parquet approach is that even the simplest approximation does 
not allow for explicit nonperturbative solutions in closed form. Hence we need further simplifications to 
obtain a workable approximate scheme. We proposed such a scheme and tested it on effective impurity 
models with frequencies as dynamical variables. The main idea of our simplification is to abandon 
all irrelevant variables that do not influence the existence and quality of singularities in two-particle 
functions. We reduced the number of relevant two-particle channels at half filling to two and kept only 
the variables at which the divergence in the appropriate channel arises. The resulting simplified theory is 
comparable in complexity with the single-channel (FLEX) approximations, but behaves in a qualitatively 
different manner at intermediate and strong coupling. We showed explicitly on the Hubbard and single- 
impurity Anderson models in an external magnetic field that all simpler approximations get unstable at 
intermediate coupling and fail to describe the transition from weak to strong coupling. Only the parquet 
approximation remained stable on the metallic side of the Mott-Hubbard metal-insulator transition. 
Although analytically stable, the simplified parquet approximation with linearly distributed integration 
mesh points became numerically unstable when used in the critical region of the Mott-Hubbard metal- 
insulator transition. This instability is of purely numerical character and does not follow from the quality 
of the parquet approximation itself. 

A picture from the parquet theory for the transition from weak to strong coupling in the mean-field 
solution of the Hubbard model has not yet been completed, but the presented analysis strongly suggests 
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in agreement with [|6|, that there cannot be a sharp, second-order Mott-Hubbard transition for finite 
couphng at zero temperature. First, we have not found a sharp transition when an external magnetic field 
is applied and the only divergence in the two-particle spin-flip function was found at the boundary to the 
fully spin-polarized solution. Second, if there were a sharp transition in the paramagnetic, spin-symmetric 
phase, then inevitably two local two-particle functions from the electron-hole and the interaction channels 
had to diverge at the transition point. It follows from our stability analysis that then a local order should 
set in. However, local two-particle functions are sums of momentum-dependent two-particle functions 
over the Brillouin zone, even in infinite dimensions. If such a sum diverges, it must be at least one 
momentum qo at which the paramagnetic solution gets unstable before the Mott-Hubbard transition is 
reached. Hence, a long-range order will conceal the Mott-Hubbard metal-insulator critical point. We 
conclude that Fermi liquid can break down only if a long-range order sets in. 

Based on the experience and understanding of the origin of singularities at intermediate and strong 
coupling, we can conclude that only theories with dynamical vertex renormalizations are capable to 
describe the transition from weak to strong coupling qualitatively correctly. The parquet approximation 
with an adequate numerical technique to iterate or diagonalize the resulting nonlinear equations offers 
a suitable analytic-numerical tool for this purpose. Especially quantum criticality is the situation where 
only the parquet approach is able to give a complete picture of the critical behavior and symmetry- 
breaking order parameters. What is to be done in particular cases is to find appropriate simplifications 
making the parquet theory viable but retaining the essential physics. 

To summarize, we have presented a general parquet construction treating one- and two-particle Green 
functions on the same footing within a single self-consistent renormalization scheme. We showed that 
instabilities are induced by divergences in two-particle functions to which symmetry-breaking order pa- 
rameters are always attributed, either normal of density type or anomalous. Further on, an instability 
in an approximate theory need not definitely indicate a transition to a long-range order. A singularity in 
one two-particle channel induces singular vertex corrections in the other channels that tend to suppress 
the long-range order. Even if a divergence appears only in one channel, as in an external magnetic field, 
a singularity in the other channels is generated and at least one another channel is needed to introduce 
a two-particle self-consistence and to keep the solution stable. The parquet approach is always to be 
used where quantum phase transitions are present and stability criteria in at least two channels may 
be broken. It is the low-temperature limit in low spatial dimensions (d = 0, 1,2) and in intermediate- 
and strong-coupling regimes in three and higher dimensions. There, the parquet approximation is the 
simplest consistent and qualitatively correct description. 

The work was supported by grant No. 202/98/1290 of the Grant Agency of the Czech Republic. 



APPENDIX 



To derive the generating thermodynamic functional for the parquet approximation we proceed as 
generally outlined in and explicitly evaluated for the parquet approximation in p^ . I. e. we introduce 
pairwise Legendre conjugate variables for the irreducible and reducible two-particle vertex functions. We 
then integrate the Bethe-Salpeter equations (0) so that 
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We use the notation introduced in Sec. ^ and particularly and write the functionals as traces over 
the active variables in each channel 
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where G*^^-* is the free two-particle propagator defined in ( |32a| ). Note that the trace in the interaction 
channel includes spin as an active variable. The logarithms with subscript denote the appropriate matrix 
multiplication used in their power series expansions. 

Integration of the other set of two-particle equation of the parquet approximation leads to a functional 
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Be aware that the operator of the bare interaction acts only between different spins, i. e. contains (So-.-o-- 
The functionals $^ are exact and only the functional is approximate, since only it depends on the 
completely irreducible vertex being replaced here by the bare interaction. 

The Luttinger-Ward generating functional is a sum of $^ and completed with second-order term 
and products of the Legendre conjugate two-particle functions A and K, to keep the variation of the 
Luttinger-Ward functional w.r.t. two-particle functions zero. We obtain 
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where we introduced in each functional its explicit dependence on the variational one- and two-particle 
functions and the bare interaction U . 

The grand potential generating the complete thermodynamics of the parquet approximation is con- 
structed from the Luttinger-Ward functional and a free-electron term with the Hartree contribution . 
We hence add also the particle densities as additional variational parameters and obtain for a fixed 

= M + f ^ 

^n[n^,ni; E, G; A,/C] = XI In [iw„ + /v - e(k) - Un^a - Ecr(k, 



-t-G^(k, iujn)^^(k, iun)^ ~ Un^ni + $[G; A, IC; U] 



(A.7) 



Here n|, n^; S, G; A, K, are independent variables and complex functions. Their physical values are cho- 
sen from stationarity of the grand potential (A.7) with respect to variations of its independent vari- 
ables/functions. There are three pairs of Legendre conjugate variational variables in ( A.7 ). The Hartree 
parameters and n^, the one-electron functions Y,„{k) and G„{k), and finally the two-particle irreducible 
and reducible functions A"^, (fc, k'\ q) and /C"^.' i^j q)j respectively. 
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